clc;clear;
theta = 2 * pi/3;
%% 可视化——等高线及三维曲面
load x.mat
load y.mat
load z.mat
% 创建网格
lamda=1852;
xx=lamda*x;
yy=lamda*y;
[X, Y] = meshgrid(xx, yy);
Z = -z;
% 绘制网面形状
figure(5)
mesh(X, Y, Z); % 使用 mesh 函数绘制网面形状
% 或者使用 surf 函数绘制带有颜色映射的网面形状
% surf(X, Y, z);

% 设置坐标轴标签和标题
xlabel('东西方向');
ylabel('南北方向');
zlabel('海水深度');

% 绘制等高线图
figure(4); % 创建新的图形窗口
contour(X, Y, z, 'Fill',"on"); % 使用 contour 函数绘制等高线图
% 指定填充颜色
cmap = flipud(parula); % 使用parula颜色映射并翻转颜色顺序
contourf(X, Y, z, 'LineColor', 'k','LevelStep', 5); % 黑色等深线
colormap(cmap); % 使用jet颜色映射

% 添加颜色条
colorbar;
% 设置坐标轴标签和标题
xlabel('东西方向');
ylabel('南北方向');
% 区域划分线
line([0, 7408], [8.000640000000000e+03, 8.000640000000000e+03], 'LineWidth', 3, 'Color', 'k');
line([5.000400000000001e+03, 5.000400000000001e+03], [0, 8.000640000000000e+03], 'LineWidth', 3, 'Color', 'k');


%% 划分区域，并作出各个区域的等高线
X1 = X(1:217, 1:136); 
Y1 = Y(1:217, 1:136);
Z1 = Z(1:217, 1:136);
% figure(1)
% mesh(X1, Y1, Z1);

X2 = X(1:217, 136:201);
Y2 = Y(1:217, 136:201);
Z2 = Z(1:217, 136:201);
figure(2)
mesh(X2, Y2, Z2);
% %指定填充颜色
% cmap = flipud(parula); % 使用parula颜色映射并翻转颜色顺序
% contourf(X2, Y2, -Z2, 'LineColor', 'k','LevelStep', 5); % 黑色等深线
% colormap(cmap); % 使用jet颜色映射
% %添加颜色条
% colorbar;
% %设置坐标轴标签和标题
% axis equal
xlabel('东西方向');
ylabel('南北方向');
zlabel('海水深度');


X3 = X(217:251, 1:201); 
Y3 = Y(217:251, 1:201);
Z3 = Z(217:251, 1:201);
figure(3)
mesh(X3, Y3, Z3);
%指定填充颜色
% cmap = flipud(parula); % 使用parula颜色映射并翻转颜色顺序
% contourf(X3, Y3, -Z3, 'LineColor', 'k','LevelStep', 5); % 黑色等深线
% colormap(cmap); % 使用jet颜色映射
% %添加颜色条
% colorbar;
% %设置坐标轴标签和标题
% axis equal
axis([0, 7408, 8.000640000000000e+03, 9260, -100, -20])  
xlabel('东西方向');
ylabel('南北方向');
zlabel('海水深度');

%% 区域一的求解
% 曲面参数
a=-16.23;b=-0.002181;c=0.0001069;d=-2.192e-06;
e=1.567e-06;f=-9.33e-07;

y0=20;z0=-16.2290;
theta=2*pi/3;
fai=pi/2-theta/2;
res=[];zans=[];
suml=0;
while y0<8000
    syms y;
    func=@(y)(-(a+c*y+f*y^2))*sin(theta/2)/sin(fai)-(y-y0)-z0/sin(fai)*sin(theta/2);
    solutions = fsolve(func,y0+100);
    y0=solutions;
    res=[res,y0];
    if length(res)>=56
        suml=suml+1.6892*y0;
    end
    z0=a+c*y0+f*y0^2;
    zans=[zans,z0];
end
suml
x0=27.7128;z0=-16.2290;y0=4000;
theta=2*pi/3;
fai=pi/2-theta/2;
res=[x0];
resy=[y0];
cutz=zans(53);
while 1
    funcx=@(x)(-(a+b*x+d*x^2)/sin(fai)-(a+b*x0+d*x0^2)/sin(fai)-(x-x0));
    solutions1 = fsolve(funcx,x0+50);
    x0=solutions1;
    res=[res,x0];
    funy=@(y)(a+b*x0+c*y+d*x0^2+e*x0*y+f*y^2-cutz);
    solutions2 = fsolve(funy,y0+10);
    y0=solutions2;
    suml=suml+y0;
    resy=[resy,y0];
    if solutions1>2000
        break;
    end
end

%% 区域二的求解
% 区域二中的平面1
p11 = [X2(1, 1), Y2(217, 1), Z2(217, 1)];
p12 = [X2(1, 1), Y2(1, 1), (Z2(1, 1)+Z2(217, 66))/2];
p13 = [X2(1, 66), Y2(217, 1), (Z2(1, 1)+Z2(217, 66))/2];
pl11 = p11 - p12;
pl12 = p11 - p13;
mian1 = cross(pl11, pl12);
alpha1 = acos(dot(mian1, [0,0,1])/sqrt(mian1(1)^2+mian1(2)^2+mian1(3)^2));
% 区域二中的平面2
p21 = [X2(1, 1), Y2(1, 1), (Z2(1, 1)+Z2(217, 66))/2];
p22 = [X2(1, 66), Y2(217, 1), (Z2(1, 1)+Z2(217, 66))/2];
p23 = [X2(1, 66), Y2(1, 1), Z2(1, 66)];
pl21 = p21 - p22;
pl22 = p21 - p23;
mian2 = cross(pl21, pl22);
alpha2 = acos(-dot(mian2, [0,0,1])/sqrt(mian2(1)^2+mian2(2)^2+mian2(3)^2));
% 画出区域二的近似平面
figure(2)
hold on
fill3([X2(1, 1), X2(1, 1), X2(1, 66)], [Y2(217, 1), Y2(1, 1), Y2(217, 1)], [Z2(217, 1), (Z2(1, 1)+Z2(217, 66))/2, (Z2(1, 1)+Z2(217, 66))/2], 'r');

figure(2)
hold on
fill3([X2(1, 1), X2(1, 66), X2(1, 66)], [Y2(1, 1), Y2(217, 1), Y2(1, 1)], [(Z2(1, 1)+Z2(217, 66))/2, (Z2(1, 1)+Z2(217, 66))/2, Z2(1, 66)], 'r');

% 等效坡长
theta_xy = atan((X2(1, 66)-X2(1, 1))/Y2(217, 1));
l = Y2(217, 1) * sin(theta_xy);
% 求解测线
% 平面1
cover = 0;
D0 = -Z2(217, 1);%最浅处的深度
angle1 = pi/2 + alpha1 - theta/2;
angle2 = pi/2 - alpha1 - theta/2;
xlast = D0/(sin(angle1)/cos(alpha1)-sin(theta/2)*tan(alpha1));
pos1 = [xlast];
while 1
    xnew = (xlast * (1/cos(alpha1) + sin(theta/2) * tan(alpha1)/sin(angle2)) + D0 * (1-cover) * sin(theta/2) * (1/sin(angle2) + 1/sin(angle1)))/(cover * tan(alpha1) * sin(theta/2) * (1/sin(angle2) + 1/sin(angle1)) - tan(alpha1) * sin(theta/2)/sin(angle1) + 1/cos(alpha1));
    if xnew > l
        break;
    end
    pos1 = [pos1, xnew];
    xlast = xnew;
end
% 平面2
cover = 0;
D0 = -(Z2(1, 1)+Z2(217, 66))/2;%最浅处的深度
angle1 = pi/2 + alpha2 - theta/2;
angle2 = pi/2 - alpha2 - theta/2;
xlast = D0/(sin(angle1)/cos(alpha2)-sin(theta/2)*tan(alpha2));
pos2 = [xlast];
while 1
    xnew = (xlast * (1/cos(alpha2) + sin(theta/2) * tan(alpha2)/sin(angle2)) + D0 * (1-cover) * sin(theta/2) * (1/sin(angle2) + 1/sin(angle1)))/(cover * tan(alpha2) * sin(theta/2) * (1/sin(angle2) + 1/sin(angle1)) - tan(alpha2) * sin(theta/2)/sin(angle1) + 1/cos(alpha2));
    if xnew > l
        break;
    end
    pos2 = [pos2, xnew];
    xlast = xnew;
end
% 求侧线总长度
l_xiebian = sqrt(Y2(217, 1)^2 + X2(1, 66)^2);
L1 = 0;
for i = 1:length(pos1)
    L1 = L1 + pos1(i) * l_xiebian / l;
end
L2 = 0;
for i = 1:length(pos2)
    L2 = L2 + (l - pos2(i)) * l_xiebian / l;
end

%% 区域三的求解
p1 = [0, Y3(1,1), Z3(1,1)];
p2 = [0,Y3(35,1), Z3(35,1)];
% 绘制区域三的近似平面
figure(3)
hold on;
fill3([0, 0, 7408, 7408], [8.000640000000000e+03, 9260, 9260, 8.000640000000000e+03], [Z3(1,1), Z3(35,1), Z3(35,1), Z3(1,1)], 'r');

l = Y3(35,1) - Y3(1,1);
h = Z3(1,1) - Z3(35,1);
alpha = atan(h/l);
cover = 0;
D0 = -Z3(1,1);
angle1 = pi/2 + alpha - theta/2;
angle2 = pi/2 - alpha - theta/2;
xlast = D0/(sin(angle1)/cos(alpha)-sin(theta/2)*tan(alpha));
pos3 = [xlast];
while 1
    xnew = (xlast * (1/cos(alpha) + sin(theta/2) * tan(alpha)/sin(angle2)) + D0 * (1-cover) * sin(theta/2) * (1/sin(angle2) + 1/sin(angle1)))/(cover * tan(alpha) * sin(theta/2) * (1/sin(angle2) + 1/sin(angle1)) - tan(alpha) * sin(theta/2)/sin(angle1) + 1/cos(alpha));
    if xnew > l
        break;
    end
    pos3 = [pos3, xnew];
    xlast = xnew;
end
L3 = length(pos3) * X2(1, 66);

%% 指标计算
% 总长度
L = L1 + L2 + L3 + suml;

%% 可视画——绘制测线
% 区域二测线
k = 1/tan(theta_xy);
pp = [5.000400000000001e+03, 8.000640000000000e+03];
b0 = -k * pp(1) + pp(2);
xxx = 5.000400000000001e+03:1:7408;
yyy = k * xxx + b0;
x_y = zeros(2,length(xxx));
x_y(1,:) = xxx;
x_y(2,:) = yyy;
delta1 = pos1 / sin(theta_xy);
figure(6)
axis([0, 7408, 0, 9260]);
hold on
temp = 0;
for i = 1:length(delta1)
    x_y_temp = x_y;
    x_y_temp(2,:) = x_y_temp(2,:) - delta1(i);
    for j =1:length(x_y_temp(1,:))
        if x_y_temp(2,j) > 8.000640000000000e+03
            temp = j;
            break
        end 
    end
    x_y_temp(:,temp:length(x_y_temp(1,:))) = [];
    plot(x_y_temp(1,:), x_y_temp(2,:), 'LineWidth', 0.7, 'Color', 'r');
end
delta2 = pos2 / sin(theta_xy);
for i = 1:length(delta2)
    plot(xxx, yyy - 8.000640000000000e+03 - delta2(i),'LineWidth', 0.7, 'Color', 'r')
end

% 区域三测线
xxx = 0:1:7408;
yyy = ones(1,length(xxx)) * 8.000640000000000e+03;
for i = 1:length(pos3)
    plot(xxx, yyy + pos3(i),'LineWidth', 0.7, 'Color', 'b')
end


% 区域一测线
listf=[];
for i=53:length(zans)
    syms x y;
    listf =[listf,a+b*x+c*y+d*x^2+e*x*y+f*y^2-zans(i)];
    % 绘制曲线
    ezplot(listf(end), [0, 5000, 0, 8000]);
    hold on;
    % 添加标题和轴标签
    xlabel('x');
    ylabel('y');
end
for i=1:length(res)
    hold on;
    plot([res(i),res(i)],[0,resy(i)],'k','LineWidth', 0.7);
end
xlabel('东西方向');
ylabel('南北方向');